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1 What are frustrated systems? 

Frustrated systems are simply systems in which the individual entities that build 
up the model (like spins, bosons, fcrmions, monomers, etc.) feel some sort of 
"frustration" in the literal sense. This means that on their search for a minimal 
energy configuration at lower and lower temperatures they are not able to satisfy 
all interactions with one another or with impurities simultaneously. 

As an example we consider a model for a directed polymer in a disordered 
environment 

H = Y^{xi-Xi+if-Y,V{xi), (1) 

i i 

^ V ' ^ V ' 

(A) (B) 

where Xi is the displacement of the i-th monomer and V{x) is a random potential. 
The first term (A), the elastic energy, tries to make the polymer straight for 
T ^ 0, the second term (B) tries to bring the monomers in favorite positions, 
for which the polymer has to bend. The monomers cannot satisfy both of these 
demands simultaneously. 

Another, more famous example arc magnetic spins (for simplicity Ising spins) 
with ferromagnetic and antiferromagnetic interactions. Consider the Hamilto- 
nian for 4 spins (e.g. an elementary plaquette of a square lattice) 

H = — (7i(T2 — (T2(J3 — Cr3(J4 + 0"4(Ti (2) 

and try to find a configuration of the Ising spins CTj = ±1 that minimizes this 
simple energy function. Naively one starts with some value for the first spin, let's 
say (Ti = +1, the first term would then imply (J2 = +1, the second (T3 = +1 and 
the third = +1. But what about the last term — here cti = cr4 + 1 is not the 
most favorable configuration. Thus it is impossible to satisfy all local interactions 
at once, this is why Toulouse[l] introduced the concept of frustration for these 
plaquette occurring naturally in spin glass models. After some thought one finds 
that many (i.e. 8) different spin configurations for (2) have a minimal energy, 
but all of them break one bond. This is the notorious frustration induced ground 
state degeneracy. 

This kind of frustration can occur either via quenched disorder (i.e. a ran- 
dom, time-independent distribution of ferromagnetic and antiferromagnetic spin 
interactions) or without any disorder, for example in the fully frustrated anti- 
ferromagnetic Ising model on a triangular lattice. Of course the same problem 
occurs with XY-spins, like the XY-antiferromagnet on a triangular or Kagome 
lattice. In this letter we treat exclusively disorder induced frustration. The de- 
termination of ground states of regularly frustrated systems usually do not need 
such algorithmic tools as discussed in this lecture (see [2] and references therein 
for a number of examples). 
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2 What is special for simulations of disordered systems? 

As we learned from our simple 4-spin Hamiltonian above, frustration is often 

responsible for the existence of many degenerate (or nearly degenerate) states 
and metastable states. Suppose one intends to perform a eonventional Monte 
Carlo simulation with single spin flip heat bath dynamics of such a system. To 
explore the whole energy landscape in order to find the most favorable configu- 
rations one has overcome large energy barriers between the various minima. As a 
consequence, the relaxation times typically become astronomically large — not 
only at a possible phase transition (in which case it would be "critical slowing 
down"), but also below and above. Thus, as is well known in the community of 
computational physicists (also among experimentalists, by the way) investigat- 
ing disordered or amorphous materials, the equilibration is nearly impossible for 
large systems. Our first commandment in this context therefore is 



To be modest in system size is mandatory! 



Of course everything would be fine if there would be an efficient^ chister 
algorithm at hand, as discussed in this school. However, there are none — with 
a few exceptional cases. To invent an efficient cluster algorithm for some model, 
one first has to have a deep understanding of its physics and a knowledge or an 
intuition about the low lying energy configurations, the excitations etc. Thus, 
cum grano salis, if you have understood the system to a rather complete extent, 
you might be ready to formulate a cluster algorithm — with which you can 
add some precise numbers (critical exponents etc.) to your basic understanding. 
Unfortunately, after several decades of research we have still not reached this 
desirable state for most of the interesting disordered systems. 

The next observation is that different samples, i.e. different disorder realiza- 
tions, can have completely different dynamical and static properties. This goes 
under the name "large sample to sample fiuctuations" , which originates in the 
lack of self averaging in some physical observables. Not all of them show this no- 
torious behavior: the ground state energy, for instance, is well behaved, simply 
because the various local minima are nearly degenerate, but e.g. spatial cor- 
relation functions or susceptibilities are not self averaging quantities. Consider 
for instance the diluted ferromagnet with site concentration c and imagine the 
following two extreme situations: On a square lattice with N spins one could 
distribute N/2 spins in such a way that 1) they form a single, compact cluster, 
or 2) they occupy a sublattice such that none of them has an occupied nearest 
neighbor site. Obviously the magnetization or susceptibility has completely dif- 
ferent characteristics in the two cases: 1) is a bulk ferromagnet of volume N/2 
and will have a tendency to order ferromagnetically at some temperature Tc (in 
the limit A'' — > oo), 2) is a collection of isolated spins that will never order. 

^ We emphasize this word, because it is easy to formulate an algorithm that constructs 
some clusters. Question is, whether the flip acceptance rate is reasonable. 
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Thus one easily recognizes that the probability distribution Pl{0) of some 
observable O is usually extremely broad, in particular non-Gaussian. Rare events 
(i.e. disorder configurations with small probability) can have a strong impact on 
averaged quantities like susceptibilities or autocorrelations. This leads to our 
second commandment for all investigations of disordered systems: 



Sample a huge (!) number of disorder configurations 



The study of the probability distribution Pl{0) can be more useful than 
only average values. 

3 What can we leeirn from ground state calculations? 

The ground state is the configuration in which the "equilibrated" system settles 
at exactly zero temperature (if there is more than one, replace state by states). 
However, T = is not accessible in the real world, so why should we bother? 
There is a number of reasons for it, some of them are listed below: 

1) As long as one is interested in equilibrium properties (and not in relaxational 
dynamics, aging, etc.) an exact ground state is more valuable than a non- 
equilibrium low temperature simulation. 

2) One might expect that some features of the ground state persist at small 
temperatures (like the domain structure in the three-dimensional random 
field Ising model [fractal or not?], etc.). 

3) If the phase transition into an ordered, may be glassy state, happens at 
T = 0, one can extract critical exponents from ground state calculations (as 
for instance in the two-dimensional Ising spin glass). 

4) If the RG (renormalization group) flow for a finite T transition is governed by 
a zero-temperature fixed point, one can again extract the critical exponents 
via ground state calculations (like in the three-dimensional random field Ising 
model). These arc then, if the RG-picture is correct, identical with those for 
the finite T transition. 

5) The zero temperature extrapolation of analytical finite T predictions for the 
glassy phase can be checked explicitly, like in the SOS (solid-on-solid) model 
with a disordered substrate. 

As a motivation this should be enough, in the next section we jump directly 
in medias res. But before we start: many people think that combinatorial opti- 
mization is essentially the traveling salesman problem, only because it is far the 
most famous problem (see [3] for an excellent introduction). This is similar to 
saying that frustrated disordered systems are essentially spin glasses (actually 
the traveling salesman problem and the spin glass problem are intimately con- 
nected via their complexity). One aim of this lecture is to remove this prejudice 
and to demonstrate that there are many more problems out there than only 



6 Heiko Rieger 



spin glasses (or traveling salesmen): algorithmically much easier to handle but 
equally fascinating. This is also the reason why on the algorithmic side we focus 
mainly on network flows: with the help of the material presented in the appendix 
everybody should be able to sit down in front of the computer and to implement 
efficiently the algorithms discussed there. If someone wants to know it all, i.e. all 
background material on graph theory, linear programming and network flows, 
we refer to standard works such as [4], [5], [6], [7], [8], [10]. 

4 Ground state interface in a random medium 

Although it was historically not the flrst random Ising model that has been 

investigated with the help of the maximum flow / minimum cut algorithm (this 
was the random field Ising model, which we shall discuss later), it might be 
pedagogically more advantageous to start with the random bond Ising model 
with a boundary induced interface. The reason for the greater intuitive appeal 
of the latter problem is that the minimum cut, which the algorithm searches, is 
identical with the minimum energy interface of the physical system, which we 
search. 

The random bond Ising ferromagnet (RBIFM) is defined by 

H = -^Jij(Ti(Tj (3) 

(y) 

with fTj = ±1 Ising spins and > ferromagnetic interactions strengths be- 
tween neighboring spins. These are random quenched variables, which means 
that they are distributed according to some probability distribution and fixed 
right from the beginning. («j) denotes nearest neighbor pairs of a d+l-dimensional 
lattice of size L'^ x H. We denote the coordinates by i = (xi, . . . , Xd, y). 

Since the interactions arc all ferromagnetic, the ground state is simply given 
by CTj = +1 for all sites i or aj = —1 for all sites i. Thus, up to now there 
is disorder, but no frustration in the problem. This changes by the boundary 
conditions (b.c.) we define now: wc apply periodic b.c. in the x-directions and 
fix the spins at y = 1 to be +1 and those at y = if to be —1. 

0-(a:i,...,x<i,y=l) = +1 and (T(xi,...,x<i,!/=H) = -1 (4) 

This induces an interface through the sample where bonds have to be broken, 
as indicated in fig. 1. If all bonds would be of the same strength = J we would 
have the pure Ising model and the interface would simply be a d-dimensional 
hyperplane perpendicular to the y-direction, which costs an energy of JL'^, J for 
each broken bond. Because of the randomness of the Jy it is energetically more 
favorable to break weak bonds: the interface becomes distorted and its shape is 
rough. This model has also been used to describe frac;tm-es in materials where 
the Jij represents the local force needed to break the material and it is assumed 
that the fracture occurs along the surface of minimum total rupture force. 
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Fig. 1. Two-dimensional random bond Ising model with a binary distribution of in- 
teraction strengths Jij € {Jj, J2} with Ji ^ J2 > 0: Thick lattice bonds are strong 
{Jij = Ji), broken lattice bonds are weak Jij — J2. The left and right boundary spins 
are connected to two extra spins as = +1 and at = —1, respectively, with infinitely 
strong bonds. \ means at = +1, \ means (7, = —1. The dotted line is the resulting 
interface: in this example it passes only weak bonds, by which it is of minimal energy. 
It partitions the lattice into up-and down-spins, the bonds (or more precisely: their cor- 
responding forward arcs — see appendix A. 3) that intersect the dotted line constitute 
the set {S,S) for the s-t-cut [3,8]. 



How do we solve the task of finding the minimal energy configuration for 
the interface? First we map it onto a flow problem in a capacitated network 
(see appendix A for the nomenclature). We introduce two extra sites, a source 
node s, which we connect to all spins of the hyperplane y = 1 with bonds 
Js,{xi,...,xd,v=i) — -^oo' ^ node t, which we connect to all spins of the hy- 
perplane y = H with bonds Js,{xu...,xa,y=H) = Joo- We choose Joo = '^J2{ij) Jio^ 
i.e. strong enough that the interface cannot pass through a bond involving one 
of the two extra sites. Now we enforce the b.c. (4) by simply fixing as = +1 and 
at = —1. The graph underlying the capacitated network we have to consider is 
now defined by the set of vertices (or nodes) 

TV = {!,..., i'^+i}u{s,0 (5) 

and the set of edges (or arcs) connecting them 

A={{i,j)\i,jGN, Jij>0}. (6) 

Note that wc have forward and backward arcs for each pair of interacting sites 
in the lattice. The capacities Uij of the arcs is given by the bond strength 
Jij. For any spin configuration a = (cti, . . . , fT]v) we define now 



S = {i€ N\ai = +1} (7) 
S={i& N\ai = -1} = N\S 
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Obviously as G S and at G S. The knowledge of S is sufficient to determine the 
energy of any spin configuration via (3): 

His) = - E _ + E _ J^^i (8) 

(i,j)£E(S) {i,j)GE(S) (i,i)e(S,S) 

= -C + 2 

where E{S) = G S,j G 5}, E(S) = {{t,3)\t G G 5} and (5,5) = 

G 5, j G S}. The constant C = J2{i j)eE{N) '^y irrelevant (i.e. indepen- 
dent of S) . Note that {S, S) is the set of edges (or arcs) connecting S with S, this 
means it cuts A'' in two disjoint sets. Since s £ S and t £ S, this is a so called 
s-f-cut-set, abbreviated [S,S]. Thus the problem of finding the ground state of 
(3) with the interface inducing b.c. (4) can be reformulated as a minimum cut 
problem 

minscw {•^'(5)} = minjg;5j ^ Jij. (9) 

{i,3)e{S,S) 

in the above defined capacitated network (with H' = [H + C)/2). It does not 
come as a surprise that this minimum cut is identical with the interface between 
the ((7j=+l)-domain and the ((7j=-l)-domain that has the lowest energy. Ac- 
tually any s-t-cut-set defines such an interface, some of them might consist of 
many components, which is of course energetically unfavorable. 

To conclude, wo have to find the minimum cut in a capacitated network, 
which is, as we show in appendix A, equivalent to finding a maximum flow 
from node s to node t. An intuitive argument for this famous max-flow-min-cut 
theorem is the following: Suppose you have to push, let's say waterflow through 
a network of pipelines, each with some c;apacity. The capacities in our case are 
the ferromagnetic interaction strengths on the bonds (pipes) between the nodes. 
Somewhere in the network there is a bottle-neck (in general consisting of several 
pipes) which does not allow a further increase of the waterflow sent from the 
source to the sink. If the maximum possible flow goes through the network, the 
flow on the pipes of the bottle-neck is equal to their capacity. The minimum cut 
is simply the global bottle-neck with the smallest capacity, and thus determines 
the maximum flow. 

One can solve the above task by applying the straightforward augmenting 
path algorithm discussed in appendix A. 2, which is based on the idea to find 
directed paths in the network on which one could possibly send more flow from 
the source to the sink. If one flnds such paths, one augments the flow along them 
(i.e. pushes more water through the pipes), if there are none, the present flow is 
optimal. In the latter case one identifies the corresponding s-t cut, which then 
yields the exact ground state interface for the above problem. 

A more efficient way is to use the preflow push algorithm presented in ap- 
pendix A. 4. The idea of this algorithm is to flood the network starting from 
the source. Then one encounters the situation that some nodes are not able to 
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transport the flood coming from the source into the direction of the sink, which 
means that one has to send some flow back to the source. The time consuming 
part of this algorithm is the retreat of floods that have been pushed too far, 
and this retreat happens faster if the capacities of backwards arcs is as large as 
possible. Bearing this observation in mind Middleton [11] has suggested a nice 
modification of the original problem that yields a significant speed up: to forbid 
overhangs of the interface we are discussing is equivalent to introduce backward 
arcs with infinite capacity in the corresponding fiow network (obviously a mini- 
mum cut will then never contain such an arc as forward arc). Thus, in the case 
that too much flow has been pushed, the retreat works with maximum efficiency. 

To conclude let us cite a number of results that have been obtained in this 
way. Of particular interest here is the width of the interface 



where yx is the j/-coordinate of the point (x, y) of the interface (note that because 
of the "no overhangs" prescription yx is single- valued). [. . .]av means an average 
over the disorder. One expects the finite size scaling form as indicated with a 
roughness exponent From the ground state calculations and finite size scaling 
one finds [11] that C = 0.41 ± 0.01 in 2d (L up to 120, and H up to 50, with 
10^ - 10^ samples), and C = 0.22 ± 0.01 in 3d {L up to 30 and H up to 20). 

5 The random field Ising model 

The random field Ising model (RFIM, for a review see [12]) is defined 



with (7i = ±1 Ising spins, ferromagnetic bonds Jij > (random or uniform), (ij) 
nearest neighbor pairs on a rf-dimensional lattice and at each site i a random field 
hi £ R that can be positive and negative. The first term prefers a ferromagnetic 
order, which means it tries to align all spins. The random field, however, tends 
to align the spins with the field which points in random directions depending on 
whether it is positive or negative. This is the source of frustration in this model. 

Let us suppose for the moment uniform interactions Jij = J and a sym- 
metric distribution of the random fields with mean zero and variance hr- It is 
established by now that in 3d (and higher dimensions) the RFIM shows ferro- 
magnetic long range order at low temperatures, provided hr is small enough. 
In Id and 2d there is no ordered phase at any finite temperature. Thus in 3d 
one has a paramagnetic/ferromagnetic phase transition along a line hc{T) in the 
/ir-T^-diagram. 

The renormalization group (RG) picture says that the nature of the transition 
is the same^ all along the line hc{T), with the exception being the pure fixed 

^ We leave aside the discussion about a possible tricritical point (which does not seem 
to be the case [12]) and the existence of an intervening spin glass phase. 



WiL,H) = {[ylU - [y.]iy^' ~ L'^wiH/L^) 



(10) 




(11) 
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point athr = and Tc ~ 4.51J. The RG flow is dominated by a zero temperature 
fixed point at hc{T=0). As a consequence, the critical exponents determining the 
critical behavior of the RFIM should be all identical along the phase transition 
line, in particular identical to those one obtains at zero temperature by varying 
hr alone. 

Therefore we consider zero temperature from now on. Close to the transition 
at he = hc{T=0) one would e.g. expect for the disconnected susceptibility 



1 

^dis — -p- 



m^-s) (12) 



where 6 = hr — he is the distance from the critical point and v is the thermal 
critical exponent. An analogous expression holds for the magnetization involv- 
ing the exponent /3. Thus to estimate a set of critical exponents the task is to 
calculate the ground state configurations of the RFIM (11). 

This optimization task is again equivalent to a maximum flow problem [13], 
as in the interface model discussed in the last section. Historically the RFIM 
was the first physical model that has been investigated with a maximum flow 
algorithm [14]. However, here the minimum-cut is not a geometric object within 
the original system and therefore we found it more intuitive to discuss the RFIM 
after the interface problem. 

In essence we proceed in the same way as in the last section. Again we add 
to extra nodes s and t and put spins with fixed values there: 

as = +1 and at = -1 (13) 

We connect all sites with positive random field to the node s and all sites with 
negative random field to t: 



J si 
Jit 



if hi >0 

ifhi<0 
if /li < 
if /li > 



(14) 



We construct a network with the set of nodes iV = {1, • • • , L'^} U {s, t} and the 
set of (forward and backward) arcs A = i,j G N, Jij > 0}. Each of them 

has a capacity Uij = Jij . Now wc can write the energy or cost function as 

E = - Jijataj (15) 

and, by denoting the set S = {i & N\Si = +1} and S = N\S the energy can be 
written as in equation (8): 

E{S) = -C + 2 J2 (16) 

(«,j)e(s,s) 
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with C = X)(j^j)gA • "^^^ problem is reduced to the problem of finding a 
minimum s-f-cut as in (9). The difference to the interface problem is that now 
the extra bonds connecting the two special nodes s and t with the original lattice 
do not have infinite capacity: they can lie in the cut, namely whenever it is more 
favorable not to break a ferromagnetic bond but to disalign a spin with its local 
random field. In the extended graph which we consider the s-f-cut again forms 
connected interface, however, in the original lattice (without the bonds leading 
to and from the extra nodes) the resulting structure is generally disconnected, a 
multicomponent interface. Each single component surrounds a connected region 
in the original lattice containing spins, which all point in the same direction. In 
other words, they form ferromagnetically ordered domains separated by domain 
walls given by the subset of the s-t-cxxt that lies in the original lattice. 

The maximum flow algorithm has been used by Ogielski [14] to calculate the 
critical exponents of the RFIM via the above mentioned finite size scaling. He 
obtained 

!y=1.0±0.1, ^=1.1 ±0.2, /3»0.05 (17) 

with j3 being so small that it is (numerically) indistinguishable from zero, in- 
dicating a discontinuous transition. These estimates are compatible with those 
obtained by recent Monte Carlo simulations supporting the RG idea of the uni- 
versality of the transition at finite and zero temperature. However, this is still 
not the end of the story, since various scaling predictions, also based on the RG 
picture, are violated. For further details we refer to the review [12]. 

6 The diluted antiferromagnet in an external field 

Experimentally it is of course hard to prepare a random field at each lattice site, 
therefore one might ask why people have been so enthusiastic about the RFIM, 
discussed in the last section, over decades. Actually it is because within a field 
theoretic perturbation theory (around small random fields) it has been shown [15] 
that the RFIM is in the same universality class as the diluted antiferromagnet 
in a uniform magnetic field (DAFF) defined via 

-ff = + ^ JijSiSjCFicrj - ^ hiSiCFi (18) 

{ij) i 

where Gi = ±1, > 0, {ij) are nearest neighbor pairs on a simple cubic lattice, 
and Si e {0, 1} with = 1 with probability c, representing the concentration 
of spins. Usually one takes Jij = J and hi = h, both uniform. Because of the 
plus sign in front of the first term in (18) all interactions are antiferromagnetic, 
the model represents a diluted antiferromagnet, for which many experimental 
realizations exist (e.g. FccZni_cF2). Now that neighboring spins tend to point 
in opposite directions due to their antiferromagnetic interaction a uniform field 
competes with this ordering tendency by trying to align them all. Thus it is 
again a frustrated system. Due to the analogy to the RFIM model one expects 
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at low temperatures and small enough fields a second order phase transition 
from a paramagnetic to an antiferromagnetic phase. 

In recent years people began to doubt the folklore that the DAFF is under 
all circumstances a good experimental realization of the RFIM model. Note that 
this result has been derived for small fields h, and the question is whether this 
still holds at larger fields. The largest field value at which the paramagnetic- 
antiferromagnetic transition can be studied is hc{T = 0). This motivates the 
study of the ground state transition along the same lines as in the RFIM context. 
Preliminary results [17] indicate that the critical exponents are different here, 
which implies that the RFIM and the DAFF are in diff'erent universality classes 
at large field values. 

Here we are primarily interested in the question whether we can again map 
the calculation of ground states onto a maximum fiow problem, as for the RFIM. 
The answer is yes as long as the antiferromagnetic interactions are short ranged 
among nearest neighbors on a bipartite lattice. With zero external field the 
ground state would be antiferromagnetic, which means (remember we have a 
simple cubic lattice) that we can define two bipartite sublattices A and B like 
the black and white fields of a checkerboard. Each site i'm A finds all its nearest 
neighbors j in B and vice versa. Define new spin and field variables via 

/ _ J +(Ji for i G A _ / for i G A , , 

\-ai iov i&B ' ' \-eihi for i&B ' ^ ' 

Since cr-tr^- = — (TjCTj for all nearest neighbor pairs (ij) one can write (18) as 

H = -Y,J[,aW,-Y,h!A (20) 

with = JijSiSj. Now the Hamiltonian has exactly the same form as the one for 
the RFIM, since J-j > 0. Note that even if one starts with uniform bonds J,j = J 
and a uniform field hi = h the dilution generates bond- and field disorder! Now 
that one has reduced the problem to the RFIM we also know how to map it to a 
maximum flow problem. Hartmann and Usadel [16] have extended the algorithm 
in such a way that all ground states can be calculated: for uniform bonds and 
fields the resulting RFIM has a discrete distribution of random bonds and fields, 
which leads in general to a high degeneracy of the ground state, something that 
does not happen in case of a uniform distribution, where usually the ground 
state is unique. 

In this context we would like to mention the Coulomb glass model [18], [19], 

which is a model for point charges on a d-dimensional lattice with long-range 
Coulomb interactions (repulsive of strength V/r with V positive and r being the 
Euclidean distance between two charges): 



(21) 
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where now the sum is over all pairs of sites of the lattice, rij € {0,1} indicates the 
presence (rij = 1) or absence (rij = 0) of a charged particle at site i and is the 
Euclidean distance between site i and site j. The random local potentials /Xj G 
[— VK, VK] represent the quenched disorder. Obviously this model is equivalent 
to an antiferromagnet with long-range interactions and random fields. Up to 
now no way of mapping this interesting problem onto a network flow problem is 
known, it seems to fall into the spin glass class, which we discuss now. 

7 The spin glass problem 

Spin glasses are the prototypes of (disordered) frustrated systems, their history 

is quite a long one and for the present status of numerical investigation I refer to 
[12], where also numerous references to experimental and theoretical introduc- 
tions can be found. In the models we discussed up to now, the frustration was 
caused by two separate terms of different physical origin (interactions and exter- 
nal fields or boundary conditions). Spin glasses are magnetic systems in which 
the magnetic moments interact ferro- or antiferromagnetically in a random way, 
as in the following Edwards- Anderson Hamiltonian for a short ranged Ising spin 
glass (SG) 

iJ = — ^ fTiCTj , (22) 
(u) 

where (7^ = ±1, (ij) are nearest neighbor interactions on a d-dimc^nsional lattice 
and the interaction strengths G R are unrestricted in sign. In analogy to eq. 
(8-9) one shows that the problem of finding the ground state is again equivalent 
to finding a minimal cut [S, S] in a network 

mmg.{H'{a)} = mm^sg^ J2 Jij, (23) 

(ij)e(s,s) 

again H' = {H + C)/2 with C = Jij- However, now the capacities (/y = Jij 

of the underlying network are not non-negative any more, therefore it is not 
a minimum-cut problem in the sense of appendix A. 3 and thus it is also not 
equivalent to a maximum flow problem, which we know how to handle eSiciently. 

It turns out that the spin glass problem is much harder than the questions 
we have discussed so far. In general (i.e. in any dimension larger than two and 
also for 2d in the presence of an external field) the problem of finding the SG 
ground state is TVP-complete [20], which means in essence that no polynomial 
algorithm for it is known (and also that chances to find one in the future are 
marginal). Nevertheless, some extremely eflicient algorithms for it have been 
developed [22], [23], [24], which have a non-polynomial bound for their worst 
case running-time but which terminate (i.e flnd the optimal solution) after a 
reasonable computing time for pretty respectable system sizes. 

First we discuss the only non-trivial case that can be solved with a polynomial 
algorithm: the two-dimensional Ising SG on a planar graph. This problem can 
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be shown to be equivalent to finding a minimum weight perfect matching, which 
can be solved in polynomial time. We do not treat matching problems and the 
algorithms to solve them in this lecture (sco [5], [6], [8]), however, we would like 
to present the idea [20] . To be concrete let us consider a square lattice with free 
boundary conditions. Given a spin configuration a (which is equivalent to —a) 
we say that an edge (or arc) (i, j) is satisfied if Jijaidj > and it is unsatisfied 
if JijCfidj < 0. Furthermore we say a plaqucttc (i.e. a unit cell of the square 
lattice) is frustrated if it is surrounded by an odd number of negative bonds (i.e. 
Jij ■ Jjk • Jki ■ Jii < with i, j, k and I the four corners of the plaquette)). 
There is a one-to-one correspondence between equivalent spin configurations (ct 
and —a) and sets of unsatisfied edges with the property that on each frustrated 
(unfrustrated) plaquette there is an odd (even) number of unsatisfied edges. See 
fig. 2 for illustration. 
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Fig. 2. Two-dimensional Ising spin glass witli ±-J couplings: Tliin lines, are positive 
interactions, thick lines are negative interactions, / means ai = +1, y means crj = — 1, 
shaded faces are frustrated plaquettes, broken lines cross unsatisfied edges. 



Note that 

if(^) = -C + 2 (24) 

unsatisfied edges 

which means that one has to minimize the sum over the weights of unsatisfied 
edges. A set of unsatisfied edges will be constituted by a set of paths (in the dual 
lattice) from one frustrated plaquette to another and a set of closed circles (see 
fig. 2). Obviously the latter always increase the energy so that we can neglect 
them. The problem of finding the ground state is therefore equivalent to finding 
the minimum possible sum of the weights of these paths between the frustrated 
plaquettes. This means that we have to match the black dots in the fig. 2 with 
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one another in an optimal way. One can map this problem on a minimum weight 
perfect matchincp problem, which can be solved in polynomial time (see [20] for 
further details). 

Note that for binary couplings, i.e. = ± J, where Jij = + J with probabil- 
ity p the weight of a matching is simply proportional to the sum of the lengths of 
the various paths connecting the centers of the frustrated plaquettes, which sim- 
plifies the actual implementation of the algorithm. In [21] the 2d ±J spin glass 
and the site disordered SG^ has been studied extensively with this algorithm. 

As we mentioned, in any other case except the planar lattice situation dis- 
cussed above the spin glass problem is TVP-hard. In what follows we would like 
to sketch the idea of an efBcient but non-polynomial algorithm [23], [25]. To 
avoid confusion with the minimum cut problem we discussed in connection with 
maximum flows one calls the problem (23) a max-cut problem (since finding the 
minimum of H is equivalent to finding the maximum of —H). 

Let us consider the vector space R^. For each cut [S, S] define x^*'"^' e R^, 

the incidence vector of the cut, by Xe ' =1 for each edge e = G {S, S) 

and Xe = otherwise. Thus there is a one-to-one correspondence between 

cuts in G and their {0, l}-incidence vectors in R^. The cut-polytope Pc{G) of G 

is the convex hull of all incidence vectors of cuts in G: Pc{G) = conv {x*-'^''^-' € 
R^ I 5 C j4}. Then the max-cut problem can be written as a linear program 

max {i^x \x€Pc{G)} (25) 

since the vertices of Pc{G) are cuts of G and vice versa. Linear programms 
usually consist of a linear cost function m^x that has to be maximized under the 
constraint of various inequalities defining a Polytope in i?" (i.e. the convex hull 
of finite subsets of i?") and can be solved for example by the simplex method, 
which proceeds from corner to corner of that polytop to find the maximum (see 
e.g. [5], [7], [8]). The crucial problem in the present case is that it is AfV-haid 
to write down all inequalities that represent the cut polytop Pc{G). 

It turns out that also partial systems are useful, and this is the essential 
idea for an efficient algorithm to solve the general spin glass problem as well as 
the traveling salesman problem or other so called mixed integer problems (i.e. 
linear programms where some of the variables x are only allowed to take on 
some integer values, like and 1 in our case) [3], [26]. One chooses a system 
of linear inequalities L whose solution set P{L) contains Pc{G) and for which 
Pc{G) = convex hull {x G P{L)\x integer}. In the present case these are < 
X <1, which is trivial, and the so called cycle inequalities, which are based on 

^ A perfect matching of a graph G = (A', A) is a set M C. A such that each node has 

only has only one edge of M adjacent to it. 
* The site disordered spin glass is defined as follows: occupy the sites of a square lattice 

randomly with A (with concentration c) and B (with concentration 1 — c) atoms. 

Now define the interactions Jij between neighboring atoms: Jij = —J if on both 

sites are A-atoms and Jij otherwise. 
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the observation that all cycles in G have to intersect a cut an even number of 
times (have a look at the cut in fig. 1 and choose as cycles for instance the paths 
around elementary plaquettes). The most remarkable feature of this set L of 
inequalities is that the separation problem^ for them can be solved in polynomial 
time: the cutting plane algorithm which, starting from some small initial set 
of inequalities, generates iteratively new inequalities until the optimal solution 
for the actual subset of inequalities is feasible. Note that one docs not solve this 
linear programm by the simplex method since the cycle inequalities are still too 
numerous for this to work efiiciently. 

Due to the insufficient knowledge of the inequalities that are necessary to 
describe Pc{G) completely, one may end up with a nonintegral solution x*. 
In this case one branches on some fractional variable Xe (i.e. a variable with 
x*f: {0, 1}), creating two subproblems in one of which x,, is set to and in the 
other Xe is set to 1. Then one applies the cutting plane algorithm recursivley 
for both subproblems, which is the origin of the name branch-and-cut. Note 
that in principle this algorithm is not restricted to any dimension, boundary 
conditions, or to the fieldless case. However, there are realizations of it that run 
fast (e.g. in 2d) and others that run slow (e.g. in 3d) and it is ongoing research 
to improve on the latter. A detailed description of the rather complex algorithm 
can be found in [26], [25]. 

The typical questions one tries to address in the context of spin glasses is: is 
there a spin glass transition at finite temperature, below which the spins freeze 
into some configuration (i.e. {(Ti)^ ^ for T < T^)- What can we do with ground 
state calculation to answer this question? Here the concept of the domain wall 
energy plays a crucial role [27]. What a finite but small temperature does is to 
destroy the ground state order by collectively flipping larger and larger clusters 
(droplets). If the energy cost for a reversal of a cluster of linear size L increases 
with L (like AE oc L** with t/ > 0) thermal fluctuation will not be able to 
destroy long range order, and thus we have a spin glass transition at finite T^. 
If it decreases (i.e. j/ < 0) long range order is unstable with respect to thermal 
fluctuations and the spin glass state occurs only at T = 0. As an example consider 
the d-dimensional pure Ising ferromagnet, for which the ground state is all spins 
up or all down. Reversing a cluster of linear size L breaks all surface bonds of 
this cluster, which means that it costs an energy AE oc L'^~^, i.e. y = d — 1 
for the pure ferromagnet. Thus the ferromagnetic state in pure Ising systems 
is stable for d > 1, which is well known. Instead of reversing spins one usually 
studies the energy difference between ground states for periodic and antiperiodic 
boundary conditions. In [28] it has been shown that 

AE ~ Ly (26) 

^ The seperation problem for a set of inequalities L consists in cither proving that a 
vector X satisfies all inequlaities of this class or to find an inequality that is violated by 
X. A linear programm can be solved in polynomial time if and only if the separation 
problem is solvable in polynomial time [9] . 
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with y = —0.281 for the 2d Ising spin glass with a uniform distribution (thus 
there is no finite T SG transition in this case). It has been speculated that in 
the ± J case for a range of concentration of ferromagnetic bonds [29] and in the 
site-random case for some concentration of A atoms [30] a spin glass phase might 
exist at non-zero temperature T > 0. This possibility has been ruled out in [21] 
with the help of ground state calculations. 

With the above mentioned branch & cut algorithm the magnetic field de- 
pendence of the ground state magnetization mi(/i) = X^i '''iOav has been 
calculated in the 2d case with a uniform coupling distribution. In [28] it has been 
shown that it obeys finite size scaling form 

mi(/i) ~ L-'^/^m{Lh}'^) (27) 

(note d = 2) with 5 = 1.481. This value is remarkable in so far as it clearly 
violates the scaling prediction S = I — y. 

Finally we would like to focus our attention on the very important concept of 
chaos in spin glasses. This notion implies an extreme sensitivity of the SG-state 
with respect to small parameter changes like temperature or field variations. 
There is a length scale A — the so called overlap length — beyond which the 
spin configurations within the same sample become completely decorrelated if 
compared for instance at two different temperatures 

Cat ■■= [{(rm+r)T{(Ti(Ji+r)T+AT]a.v ~ exp(^-r / X{AT)^ . (28) 

This should also hold for the ground states if one slightly varies the interaction 
strengths Jij in a random manner with amplitude 6. Let a be the ground state 
of a sample with couplings Jij the ground state of a sample with couplings 
Jij + SKij , where the Kij are random (with zero mean and variance one) and d 
is a small amplitude. Now define the overlap correlation function as 

Cs{r) := [aiai+ra'ia[+,U ~ c(rjVf) , (29) 

where the last relation indicates the scaling behavior we would expect (the over- 
lap length being A ~ (5~^/^) and C is the chaos exponent. In [28] this scaling 
prediction was confirmed with 1/(^ = 1.2±0.1. 

8 The SOS-model on a disordered substrate 

Up to now we have considered Ising models exclusively. Quite recently it has been 
shown [31], [32] that many more frustrated systems are amenable to ground 
state studies of the kind we discussed so far. Consider a solid-on-solid model 
with random offsets, modeling a crystalline surface on a disordered substrate as 
indicated in fig. 3. It is defined by the following Hamiltonian (or energy function): 



H = Y,f{hi-hj) 



(30) 
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where {ij) are nearest neighbor pairs on a d-dimensional lattice {d = 1,2) and 
f{x) is an arbitrary convex {f"{x) > 0) and symmetric (/(x) = f{—x)) function, 
for instance f{x) = x^. Each height variable hi = di + rii is the sum of an integer 
particle number which can also be negative, and a substrate offset c?j € [0, 1[. For 
a flat substrate, = for all sites i, we have the well known SOS-model [34]. 
The disordered substrate is modeled by random offsets di € [0, 1[ [33], which are 
distributed independently. 




Fig. 3. The SOS model on a disordered substrate. The substrate heights are denoted 
by di e [0, 1], the number of particle on site i by £ Z, which means that they could 
also be negative, and the total height on site i hy hi = di + rii 



The model (30) has a phase transition at a temperature Tc from a (thermally) 
rough phase for T > Tc to a super-rough low temperature phase for T < Tg. 
In two dimension "rough" means that the height-height correlation function 
diverges logarithmically with the distance (7(r) = [{{hi — hi+r)^)]av = a-T-log(r) 
(with a = I/tt for f{x) = x^), "super-rough" means that either the prefactor 
on front of the logarithm is significantly larger than o • T, or that C(r) diverges 
faster than log(r), e.g. C(r) oc log^(r). 

A part of the motivation to study this model thus comes from its relation to 
flux lines in disordered superconductors, in particular high-Tc superconductors: 
The phase transition occurring for (30) is in the same universality class as a 
flux line array with point disorder defined via the two-dimensional Sine-Gordon 
model with random phase shifts 

{ij) i 

where € [0, 27r[ are phase variables, 9i G [0, 27r[ are quenched random phase 
shifts and A is a coupling constant. One might anticipate that both models (30) 
and (31) are closely related by realizing that both have the same symmetries (the 
energy is invariant under the replacement rii ^ rii + m {(pi ^ <j)i + 2iTm) with m 
an integer). Close to the transition one can show that all higher order harmonics 
apart from the one present in the Sine-Gordon model (31) are irrelevant in a 
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field theory for (30), which establishes the identity of the universality classes^. 

To calculate the ground states of the SOS model on a disordered substrate 
with general interaction function f{x) we map it onto a minimum cost flow 
model. Let us remark, however, that the special case f{x) = \x\ can be mapped 
onto the interface problem in the random bond Ising ferromagnet in 3d with 
coZumnar disorder [35] (i.e. all bonds in a particular direction are identical), by 
which it can be treated with the maximum flow algorithm wc know already. 

We deflne a network G by the set of nodes N being the sites of the dual 
lattice of our original problem and the set of directed arcs A connecting nearest 
neighbor sites (in the dual lattice) and If we have a set of height 

variables rii we deflne a flow x in the following way: Suppose two neighboring 
sites i and j have a positive (!) height difference ni — nj>0. Then we assign the 
flow value Xij = rii — rij to the directed arc [i,]) in the dual lattice, for which 
the site i with the larger height value is on the right hand side, and assign zero 
to the opposite arc (j, i), i.e. Xji = 0. And also Xij = whenever site i and j are 
of the same height. See flg. 4 for a visualization of this scheme. Obviously then 
we have: 

Vi G TV : = 2;,. . (32) 

{j|W)eA} {j|(j,i)eA} 

On the other hand, for an arbitrary set of values for Xij the constraint (32) has to 
be fulflUed in order to be a flow, i.e. in order to allow a reconstruction of height 
variables out from the height differences. This observation becomes immediately 
clear by looking at flg. 4. 

We can rewrite the energy function as 

il(x) = Y hij{xij) , with Cij{x) = f{x - dij) , (33) 

with dij = di — dj. Thus our task is to minimize i?(x) under the constraint (32), 
which is (see appendix C.5) a minimuTn-cost-flowpvohlem with the mass balance 
constraints (32) and arc convex cost functions hij (xij ) . It is worth mentioning 
that this mapping from the SOS model to a flow problem is closely related to 
the dual link representation of the XY-model in two dimensions [36] . This does 
not come as a great surprise since wc already pointed out the relationship with 
the Sine-Gordon Hamiltonian involving phase variables (31). 

The most straightforward way to solve this problem is to start with all height 
variables set to zero (i.e. x = 0) and then to look for regions (or clusters) that can 
be increased collectively by one unit with a gain in energy. This is essentially 
what the negative cycle canceling algorithm discussed in appendix C.5 does: 
The negative cycles in the dual lattice surround the regions in which the height 
variables should be increased or decreased by one. However, it turns out that this 
is a non-polynomial algorithm, the so called successive shortest path algorithm is 
more efficient and solves this problem in polynomial time, see appendix C.5. This 

® Note, however, that that far away from Tc, as for instance at zero temperature, there 
might be differences in the two models. 
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Fig. 4. The flow representation of a surface (here a "mountain" of height rii — 3). 
The broken hncs represent the original lattice, the open dots are the nodes of the dual 
lattice. The arrows indicate a flow on the dual lattice, that results from the height 
differences of the variables rii on the original lattice. Thin arrows indicate a height 
difference of Xij = 1, medium arrows Xij = 2 and thick arrows Xij = 3. According to 
our convention the larger height values are always on the right of an arrow. Observe 
that on each node one the mass balance constraint (32) is fulfilled. 

algorithm starts with an optimal solution for H{x), which is given by Xij = +1 
for dij > 1/2, Xij = —1 for dij < —1/2 and Xij = for dij € [—1/2, +1/2]. 
Since this set of flow variables violates the mass balance constraints (32) (in 
general there is some imbalance at the nodes) the algorithm iteratively removes 
the excess/deficit at the nodes by augmenting flow. 

Let us briefly summarize the results one obtains by applying this algorithm 
to the ground state problem for the SOS model on a disordered substrate [32]: 

— The height- height correlation function diverges like C(r) oc log^(r) with the 
distance r. 

— XL = - f^j)\v can nicely be fltted to Xi = a + 61og(L) + 
clog^(L), indicating again a log^ dependence of the height-height correlation 
function. Moreover, the coefiicients a, b and c depend on the power n in 
f{x) = c increases systematically with increasing n. 

— By considering a boundary induced step in the ground state configuration 
one sees that the step energy increases algebraically with the system size: 
AE ~ with y = 0.45 ± 0.05. This corresponds to the domain wall energy 
introduced in the context of spin glasses in the las section. Furthermore the 
step is fractal with a fractal dimension close to 3/2. 

— Upon a small, random variation of the substrate heights di of amplitude S 
the ground state configuration decorrelates beyond a length scale X ^ 
with T] = 0.95 ± 0.05. This implies the chaotic nature of the glassy phase in 
this model in analogy to spin glasses. 

We would like to mention that this mapping of the original SOS model (30) 



Frustrated Systems: Ground State Properties 



21 



on the flow problem works only for a planar graph (i.e. free or fixed boundary 
conditions), otherwise it is not always possible to reconstruct the height variables 
m from the height differences Xij. As a counterexample in a toroidal topology 
(periodic boundary conditions) consider a flow, which is zero everywhere except 
on a circle looping the torus, where it is one. Although this flow fulfills the mass 
balance constraints (which are local) it is globally inadmissible: To the right of 
this circle the heights should be one unit larger than on left, but left and right 
become interchanged by looping the torus in the perpendicular direction, which 
causes a contradiction. If one insists on periodic boundary conditions, which 
have some advantages due to the translational invariance, one should recur to 
the special case f{x) = \x\, which can be treated differently, as we mentioned in 
the beginning of this section. 

9 Vortex glasses and trcifRc 

Finally we would like to focus on some further applications of the minimum 
cost flow algorithms that we discussed in the last section. Since we deal with 
network flow problems it should not come as a surprise, that a number of physical 
problems involving magnetic flux lines can be mapped onto them. We already 
mentioned the Sine-Gordon model with random phase shifts (31) describing a 
flux line array with point disorder and which is related to the SOS model on a 
disordered substrate. This relationship can be made more concrete with the help 
of the triangular Ising SOS model as discussed in [35]. 

The gauge glass model describes the vortex glass transition in three-dimensional 
superconductors. If one includes the screening of the interactions between vor- 
tices one can show that in the strong screening limit, the model Hamiltonian (in 
the link representation) acquires the form [38] 



where Xij are integer vortex variables living on the links (i, j) of the dual of 
the original simple cubic lattice. They represent magnetic flux lines, by which 
they have to be divergenceless — which means that they have to fulflll the mass 
balance constraint (32). The quenched random variables bij also fulfill the same 
constraint (they have to be constructed as a lattice curl from a quenched vector 
potential). Moreover one has periodic boundary condition. 

It has been shown that this model has a vortex glass transition at zero 
temperature. Thus, for the characterization of the critical behavior either low 
temperature Monte Carlo simulations [38] or ground state calculations become 
mandatory. The latter program has been performed in a tentative way in [39] 
with a stochastic, non-exact method for small system sizes (L < 4) . The problem 
of minimizing (34) under the above mentioned constraints is a convex cost flow 
problem that can be solved in a straightforward manner with the algorithms 
presented in appendix C.5. Work in this direction is in progress [37]. 




(34) 
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A further application of the minimum cost flow algorithms with convex cost 
functions is traffic flow, which became a major research topic in physics quite re- 
cently [40] . Network flow problems naturally occur in any transportation system: 
what is the shortest path between point A and point B in a road network (short- 
est path problem), how many vessels does a steamship company need to have 
in order to deliver perishable goods between several different origin-destination 
pairs (maximum flow problem) or what is the flow that satisfies the demands at 
a number of warehouses from the available supplies at a number of plants and 
that minimizes its shipping cost (typical transportation problem = minimum 
cost flow problem). 

All of the above problems are linear problems. Whenever system congestion 
or queuing effects have to be taken into account in the model describing a "real" 
network flow, the introduction of nonlinc;ar c;osts (since queuing delays vary 
nonlinearily with flows) are mandatory. In road networks, as more vehicles use 
any road segment, the road becomes increasingly congested and so the delay on 
that road increases. For example, the delay on a particular road segment, as a 
function of the flow x on that road, might be ax/{u — x). In this expression u 
denotes a theoretical capacity of the road and a is another constant: As the flow 
increases, so docs the delay; moreover, as the flow x approaches the theoretical 
capacity of that road segment, the delay on the link becomes arbitrarily large. 
In many instances, as in this example, the delay function on each road segment 
is a convex fimction of the road segment's flow, so finding the flow plan that 
achieves the minimum overall delay, summed over all road segments, is a convex 
cost network flow model. 

It should have become clear at the end of this lecture that frustrated, disor- 
dered systems and network flows are strongly related, or even completely equiva- 
lent. The quenched disorder occurring in the physical models we discussed so far 
find their counterpart in arc capacities and costs in flow problems. Thus, many 
"daily life" networks, like transportation systems or urban traffic flows, share 
many features with disordered or even glassy systems. For instance the concept 
of chaos we encountered in spin glasses as well as in the random solid-on-solid 
model should also be valid in traffic networks: the slightest random (i.e. uncon- 
trollable) change in the capacities of the roads, as for instance after a heavy 
rain or snowfall, or locally by several accidents, should completely change the 
traffic flow pattern beyond a particular length scale. A systematic study of these 
issues is certainly of great interest, not only for the theoretical understanding of 
the intrinsically chaotic nature of complex network flows but also for practical 
reasons. 
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Concepts in network flows and basic algorithms 

In this appendix we introduce the basic definitions in the theory of network 
flows, which arc needed in the main text. It represents a very condensed version 
of some chapters of the exceUent book Network flows by R. K. Ahuja, T. L. 
Magnati and J. B. OrUn[10]. The content of the subsequent chapters is self- 
contained, so that it should be possible for the reader to understand the basic 
ideas of the various algorithms. In principle he should even be able to devise a 
particular implementation of one or the other code, although I recommend to 
consult existing public-domain (!) software libraries (e.g. [41]) first. 



A The maximum flow / minimum cut problem 
A.l BeisIc definitions 

A capacitated network is a graph G = {N, A) , where N is the set of nodes and 
A the set of arcs, with nonnegative capacities Uij (which can also be infinite) 
associated with each arc G A. In our first example of the random bond 
Ising model N is simply the set of lattice sites (plus two extra nodes, sec fig. 
1), A the bonds between interacting sites and Uij the ferromagnetic interaction 
strengths. Note that Uij > is essential. Let n = #-A/' be the number of nodes 
in G and m = =ffA the number of arcs. 

The arc adjacency list is the set of arcs emanating from a node: A{i) = 
{{i,j)\{i,j)€A}. 

One distinguishes two special nodes of N: the source node s and the sink node 
t. 

A flow in the network G is a set of nonnegative numbers Xij (or a map x : ^ — > 
R+ U {0, oo}) subject to a capacity constraint for each arc 

0<Xij<Uij y{i,j)eA. (35) 

and to a mass balance constraint for each node 




s 

t (36) 



This means that at each node everything that goes in has to go out, too, with 
the only exception being the source and the sink. What actually flows from s to 
t'lsv, the value of the flow. 

The maximum flow problem for the capacitated network G is simply to 
find the flow x that has the maximum value v under the constraint (35) and 
(36). 

We make a few assumptions: 1) the network is directed, which means that 
for instance is an arc pointing from node i to node j , 2) whenever an 
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arc belongs to a network, the arc also belongs to it or is added with 
zero capacity, 3) all capacities are nonnegative integers, 4) the network does not 
contain a directed path from node s to node t composed only of infinite capacity 
arcs, 5) the network does not contain parallel arcs7 

A. 2 Residual Network and generic augmenting path algorithm 

Now that we have defined the maximum flow problem, we have to introduce 

some tools with which it can be solved. The most important one is the notion 
of a residual network, which, as it is very often in mathematics, is already half 
the solution. If we have found a set of numbers x that fulfill the mass balance 
constraints, we would like to know whether this is already optimal, or on which 
arcs of the network we can improve (or augment in the jargon of combinatorial 
optimization) the flow. 

Given a flow x, the residual capacity of any arc £ A is the maximum 
additional flow that can be sent from node i to node j using the arcs and 
(j, i). The residual capacity has two components: 1) Uij—Xij, the unused capacity 
of arc 2) xji the current flow on arc {j, i), which we can cancel to increase 

the flow from node i to j. 

The residual network G(x) with respect to the flow x consists of the arcs with 
positive residual capacities. 

An augmenting path is a directed path from the node s to the node t in the 
residual network. The capacity of an augmenting path is the minimum residual 
capacity of any arc in this path. 

Obviously, whenever there is an augmenting path in the residual network 
G(x) the flow X is not optimal. This motivates the following generic augmenting 
path algorithm. 

algorithm augmenting path 
begin 
x:=0 

while G{x) contains a directed path from node s to t do 
begin 

identify an augmenting path P from node s to node t 

5 = mm{rij\{i,j) G P} 

augment 5 units of flow along P and update G(x) 
end 
end 

All of these assumptions can be fulfilled in the physical problems we consider by 
appropriate modifications. E.g. number 3) can be fulfilled by rescaling the bond 
strengths Jij with a factor and chopping off the decimal digits. 
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Further below we will see that the flow is indeed maximal if there is no augment- 
ing path left. The main task in an actual implementation of this algorithm would 
be the identification of the directed paths from s to t in the residual network. 

Before wc come to this point we have to make the connection to the minimum 
cut problem that is relevant for the physical problems discussed in the main 
text. 

A. 3 Cuts, labeling algorithm and max-flow-min-cut theorem 

A cut is a partition of the node set N into two subsets S and S = N\S denoted 
by [S, S]. We refer to a cut as a s-t-cut if s G 5 and t <E S. 
The forward arcs of the cut [S, S] are those arcs G A with i £ S and j G S, 
the backward arcs those with j £ S and i G S. The set of all forward arcs of 
[5,5] is denoted (5,5). 

The capacity of an s-t-cut is defined to be u[S, 5] = J2(i j)eis s) ^^o- Note that 
the sum is only over forward arcs of the cut. 

The minimum cut is a s-t-cut whose capacity u is minimal among all s-t-cuts. 

Let X be a flow, v its value and [5, 5] an s-t-cut. Then, by adding the mass 
balances for all nodes in 5 we have 

i; = ^| Xij- J2 ^3i\= Yl . (38) 

'es ^{j\{i,j)eA(i)} {j\U,i)eA{i)} (i,j)e(s,s) {i,j)e(s,s) 

Since Xij < Uij and xji > the following inequality holds 

v< Uij=u[S,S] (39) 

(i,i)e(s,s) 

Thus the value of any flow x is less or equal to the capacity of any cut in the 
network. If we discover a flow x whose value equals to the capacity of some cut 
[5, 5], then x is a maximum flow and the cut is a minimum cut. The following 
implementation of the augmenting path algorithm constructs a flow whose value 
is equal to the capacity of a s-i-cut it defines simultaneously. Thus it will solve 
the maximum flow problem (and, of course, the minimum cut problem). 

As we have mentioned, our task is to flnd augmenting paths in the residual 
network. The following labeling algorithm uses a search technique to identify a 
directed path in G(x) from the source to the sink. The algorithm fans out from 
the source node to find all nodes that are reachable from the source along a 
directed path in the residual network. At any step the algorithm has partitioned 
the nodes in the network into two groups: labeled and unlabeled. The former are 
those that the algorithm was able to reach by a directed path from the source, 
the latter are those that have not been reached yet. If the sink becomes labeled 
the algorithm sends fiow along a path (identified by a predecessor list) from s 
to t. If all labeled nodes have been scanned and it was not possible to reach the 
sink, the algorithm terminates. 
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algorithm labeling 
begin 

label node t 

while node t is labeled do 
begin 

unlabel all nodes 

set pred{j) = f or each j G N 

label node s and set list := {s} 

while list ^ and node t is unlabeled do 

begin 

remove a node i from list 

for each arc {i,j) G A{i) in the residual network do 
if rij > and node j is unlabeled then 
set pred{j) = i 
label node j 
add node j to list 

end 

if node t is labeled then augment 
end 
end 

procedure augment 
begin 

Use the predecessor labels to trace back from the sink to 
the source to obtain an augmenting path P from s to i 
S = min {rij\{i,j) e P} 

augment S units of flow along P, update residual capacities 
end 

Note that in each iteration the algorithm either performs an augmentation or 
terminates because it cannot label the sink. In the latter case the current flow is 
a maximum flow. To sec this, suppose that at this stage S is the set of labeled 
nodes and S = N\S is the set of unlabeled nodes. Clearly s G 5 and t G S. Since 
the algorithm cannot label any node in 5 from any node in S, r^j = for each 
G {S,S), which implies with (37) — Uij + Xji. Thus Xij = Uij (since 
< Xij < Uij) for all {i,j) G {S, S) and Xij = for all {i,j) G {S, S). Hence 

v= ^ Xij - ^ Xij = ^ Xij = u[S, S] . (40) 
(z,j)e{s,s) {z,j)e(s,s) {i,j)e{s,s) 

This means that the flow x equals the capacity of the cut [S, 5] , and therefore 
X is a maximum flow and [S, S] is a minimum cut. 

From these observation let us note the following conclusions: 
Max-flow-min-cut theorem: The maximum value of the flow from a source 
node s to a sink node f in a capacitated network equals the minimum capacity 
among all s-t-cuts. 
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Augmenting path theorem: A flow x* is a maximum flow if and only if the 
residual network G(x*) contains no augmenting path. 

Integrality theorem: If all arc capacities are integer, the maximum flow prob- 
lem has an integer maximum flow. 

Let n be the number of nodes, m the number of arcs and U = max{uy}. 
Since any arc is at most examined once and the cut capacity is at most nU 
the complexity of this algorithm is 0{nmU) (note that the flow increases at 
least by 1 in each augmentation). Because of the appearance of the number U 
it is a pseudo-polynomial algorithm. The so called preflow-push algorithms we 
discuss now are much more efficient, in particular they avoid the delay caused 
notoriously by some bottleneck situations. 

A. 4 Preflow-push algorithm 

The inherent drawback of the augmenting path algorithms is the computation- 
ally expensive operation of sending flow along a path, which requires 0{n) time 
in the worst case. The preflow-push algorithms push flow on individual arcs in- 
stead of augmenting paths. They do not satisfy the mass balance constraints at 
intermediate stages. This is a very general concept in combinatorial optimiza- 
tion: algorithms either can operate within the space of admissible solutions and 
try to flnd optimality during iteration, or they can fulflU some sort of optimality 
criterion all the time and strive for feasibility. Augmenting path algorithms are 
of the first kind, preflow-push algorithms of the second. The basic idea is to flood 
the network from the source pushing as much flow as the arc capacities allow 
into the network towards the sink and then reduce it successively until the mass 
balance constraints are fulfilled. 

A preflow is a function :x. : A ^ R that satisfies the flow bound constraint 
Xij < Uij and the following relaxation for the excess e{i) of each node i: 



It is e{t) > and only e(s) < 0. One denotes a node i to be active if its excess 
is strictly positive e(i) > 0. 

As mentioned, preflow-push algorithms strive to achieve feasibility. Active 
nodes indicate that the solution is infeasible. Thus the basic operation of the 
algorithm is to select an active node and try to remove its excess by pushing 
flow to its neighbors. Since ultimately we want to send flow to the sink, we push 
flow to the nodes that are closer to the sink. This necessitates the use of distance 
labels: 

We say that a distance function d : N ^ Z+ U {0} is valid with respect to a 
flow X, if it satisfies 

a) d{t) = and 

b) d{i) < d{j) + 1 for every arc (i, j) in the residual network G(x). 



01(i,i)e^} 



^ Xij > Vz G N\{s, t} . 
01(iJ)eA} 



(41) 
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If d{-) is valid then it has also the following properties (where n is the number 
of nodes): 

1) d{i) < length of the shortest directed path from node i to t in G(x) 

2) d{s) > n => G(x) contains no directed path from s to t. 
Furthermore we say that d{-) is exact if in 1) the equality holds. 
Finally an arc is admissible if d{i) = d{j) + 1. 

In the preflow-push algorithm wc push flow on these admissible arcs. If the 
active node that we are currently considering has no admissible arcs, we increase 
its distance label so that we create at least one admissible arc. 

algorithm preflow-push 
begin 

preprocess 

while the network contains an active node do 
begin 

select an active node i 
push/relabel(i) 
end 
end 

procedure preprocess 
begin 

X := 

compute the exact distance labels d{i) (1) 
Xsj = Usj for each arc {s,j) € A 
d{s) = n 
end 

procedure push/relabel (i) 
begin 

if the network contains an admissible arc {i.j) then 

push 6 = mm {e{i) , rij} units of flow from node i to j 
else 

replace d{i) by mm {d{j) + l|(i,i) € A and > 0} 

end 

Ad (1): To compute the exact distance labels we have to calculate the short- 
est distances from node t to every other node, which we learn how to do in the 
next section. 

The algorithm terminates when the excess resides at the source or at the sink, 
implying that the current preflow is a flow. Since d{s) = n after peprocessing, 
and d{i) is never decreased in push/relabel(i) for any i, the residual network 
contains no path from s to t, which means according to 2) above that there is 
no augmenting path. Thus the flow is maximal. 

As in the context of the max-flow-min-cut theorem of the last section it might 
also here be instructive to visualize the generic preflow-push algorithm in terms 
of a network of (now flexible) water pipes representing the arcs with joints being 
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the nodes. The distance function, which is so essential in this algorithm, mea- 
sures how far nodes are above the ground. In this network we wish to send water 
from the source to the sink. In addition we visualize flow in admissible arcs as 
water flowing downhill. Initially, wc move the source node upward, and water 
flows to its neighbors. In general, water flows downwards to the sink; however, 
occasionally flow becomes trapped locally at a node that has no downhill neigh- 
bors. At this point we move the node upward, and again water flows downhill 
to the sink. Eventually, no more flow can reach the sink. As we continue to 
move nodes upward, the remaining excess flow eventually flows back towards 
the source. The algorithm terminates when all the water flows either into the 
sink or flows back to the source. 

The complexity of this algorithm turns out to be O(n^m), the so called FIFO 
preflow-push algorithm, which we do not discuss here, has a complexity of 0{n^). 

B Shortest path problems 
B.l Dijkstra's algorithm 

Given a network G = (TV, A) and for each arc (i, j) € A a non-negative arc-length 
dj . In the above problem, where we had to find the exact distance labels in the 
preflow-push algorithm it is simply Cy = 1 for all arcs in the residual network. 

The task is to find the shortest paths from one particular node s to all other 
nodes in the network. Dijkstra 's algorithm is a typical label-setting algorithm to 
solve this problem (with complexity 0{n^). It provides distance labels d{i) with 
each node. Each of these is either temporarily (defining a set S) or permanently 
(defining a set S = N\S) labeled during the iteration, and as soon as a node is 
permanently labeled, d{i) is the shortest distance. The path itself is reconstructed 
via predecessor indices. 

First note that d{j) = d{i) + Cij for each arc in a shortest path from 
node s to node i, and that d{j) > d{i) + Cij otherwise. By fanning out from node 
s the algorithm uses this criterion to find successively the shortest paths. 

algorithm Dijkstra 
begin 

5:=0, S = N 

d{i) := oo for each node i €. N 
d{s) := and pred{s) := 
while |5| < n do 
begin 

let i £ S he a node for which d{i) = m.\n.{d{j)\j € S} 
5:=5U{i}, S:=S\{i} 
for each (i,j)s^(z) do 
if d{j) > d{i) + Cij then 

d{j) := d{i) + Cij and pred{j) := i 

end 
end 
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The fact that we always add a node i € S with minimal distance label d{i) en- 
sures that d{i) is indeed a shortest distance (there might be other shortest paths, 
but none with a strictly shorter distance). There are special implementations of 
this algorithm that have a much better running time than O(n^). 

B.2 Label correcting algorithm 

As we said, Dijkstra's algorithm is a label-setting algorithm. The above men- 
tioned criterion 

d{i) shortest path distances <^ d{j) < d{i) + Cij € A 

gives also rise to a so called label- correcting algorithm. 
Let us define reduced arc length (or reduced costs) via 

cf^:=c,,+d{z)-d{j). (42) 

As long as one reduced arc lengths is negative, the distance labels d{i) are not 
shortest path distances: 

d(-) shortest path distances ^ ^fj^^ V(i,j)eA (43) 

For later reference we also note the following observation. For any directed cycle 
W one has 

E 4= E (44) 

{i,j)ew {i,j)ew 

The criterion (43) suggests the following algorithm for the shortest path problem: 

algorithm label-correcting 
begin 

d{s) and pred{s) := 

d{j) := 00 for each node j £ A''\{s} 

while some arc {i,j) satisfies d{j) > d{i) + Cij icfj<0) do 
begin 

dU) := d{i) + Cij cfj = 0) 
pred{j) = i 
end 
end 

The generic implementation of this algorithm has a running time 0(min{n^mC, 
m2"}) with C = max I Cij I, which is pseudo-polynomial. A FIFO implementation 
has complexity 0{nm). 

This algorithm also works for the cases in which some costs Cij are nega- 
tive, provided there are no negative cycles, i.e. closed directed paths W with 
j)£W ^ij ^- '^^^^ instruction d{j) := d{i) + Cij would decrease 

some distance labels ad (negative) infinitum. 

If there are negative cycles, one can detect them with an appropriate mod- 
ification of the above code: One can terminate if d{k) < —nC for some node k 
(again C = max |cy |) and obtain these negative cycles by tracing them through 
the predecessor indices starting at node k. This will be useful in the next section. 
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C Minimum cost flow problems 
C.l Definition 

Tlie next flow problem we discuss combines features of the maximum-flow and 
the shortest paths problem. The algorithm that solves it therefore also makes 
use of the ideas we presented so far. Let C = {N,A) be a directed network 
with a cost Cij and a capacity Uij associated with every arc (z, j) G A. Moreover 
we associate with each node i G N a number b{i) which indicates its supply or 
demand depending on whether b{i) > or b{i) < 0. The minimum cost flow 
problem is 

Minimize z{x) = CijXij (45) 

(.i,j)eA 

subject to the mass balance constraints 

E ^^1- E ^.^=Ki) Vie AT (46) 

{j\ih3)eA} {j\ij,i)eA} 

and the capacity constraints 

0<Xij<Uij y{i,j)€A (47) 

Again wc make a few assumptions: 1) All data (cost, supply/demand, capac- 
ity) are integral^, 2) the network is directed, 3) J2i H^) t;he minimum cost 
flow problem has a feasible solution (that means, one can find a flow Xij that 
fulfills the mass balance and capacity constraints^, 4) it exists an uncapacitatcd 
directed path between every pair of nodes, 5) all arc costs are non-negative 
(otherwise one could appropriately define a revered arc). 

Again the residual network G'(x) corresponding to a flow x will play an 
essential role. It is defined in the same way as in the maximum flow problem, 
in addition the costs for the backwards arcs are reversed: a fiow Xij on arc 
£ A with capacity ?iy and cost will give rise to the arcs and (j.i) 
with residual capacities r,j = Uij — Xij and rji = Xij , respectively and costs Cij 
and —Cij respectively. 

C.2 Negative cycle canceling algorithm 

First wc formulate a very important intuitive optimality criterion, the negative 
cycle optimality criterion: A feasible solution x* is an optimal solution of the 
minimum cost flow problem, if and only if the residual network G(x*) contains 
no negative cost cycle. 

The proof is easy: Suppose the flow x is feasible and G'(x) contains a negative 
cycle. The a flow augmentation along this cycle improves the function value z{x), 

® Here the same remark holds as for the maximum flow problem, previous footnote. 
® In the physical models we discuss it is b{i) = anyway, implying x = as a feasible 
solution. 
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thus X is not optimal. Now suppose that x* is feasible and G(x*) contains no 
negative cycles and let x" ^ x* be an optimal solution. Now decompose x° — x* 
into augmenting cycles, the sum of the costs along these cycles is c • x° — c • x*. 
Since G'(x*) contains no negative cycles c ■ x" — c • x* > 0, and thcrforc c ■ x° = 
c • X* because optimality of x* implies c • x° < c • x* . Thus x° is also optimal. 

The following algorithm iterativly cancels negative cycles until the optimal 
solution is reached. 

algorithm cycle canceling 
begin 

establish a feasible flow x (1) 

while G{x) contains a negative cycle do 

begin 

use some algorithm to identify a negative cycle W (2) 

S := min {rij\{i, j) G W} 

augment 6 units of flow in the cycle cind update G(x) 
end 
end 

Ad (1): Although, as wc mentioned, in the physical problems wc discuss a 
feasible solution is obvious in most cases (e.g. x = 0) we note that in principle 
one has to solve a maximum flow problem here: One introduces two extra-nodes 
s and t (source and sink, of course) and 
Vi : b{i) > add a source arc (s, i) with capacity Usi = b{i) 
Mi : b{i) < add a sink arc {i, t) with capacity uu = —b{i). 

If the maximum flow from s to t saturates all soiucc arcs (remember b{i) = 0) 
the minimum cost flow problem is feasible and the maximum flow x is a feasible 
flow. 

Ad (2): For negative cycle detection in the residual network G{x) one can 
use the label-correcting algorithm for the shortest path problem presented in the 
last section. 

The running time of this algorithm is 0{mCU), where C maxjcy l and 
U = maxUjj, which means that it is pseudopolynomial. In the next section we 
present an alternative and more efflcient way to solve the minimum cost flow 
problem. 

C.3 Reduced cost optimality 

Remember that when wc considered the shortest path problem we introduced 
the reduced costs and obtained the shortest path optimality condition cfj = 
Cij + d{i) — d{j) > 0. This means 

— cfj is an optimal "reduced cost" for arc in the sense that it measures 
the cost of this arc relative to the shortest path distances. 

— With respect to the optimal distances, every arc has a nonnegative cost. 

— Shortest paths use only zero reduced cost arcs. 
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— Once we know the shortest distances, the shortest path problem is easy to 
solve: Simply find a path from node s to every other node using only zero 
reduced cost arcs. 

The natural question arises, whether there is a similar set of conditions for more 
general min-cost flow problems. The answer is yes as we show in the following. 
For the network defined in the last section associate a real number unre- 
stricted in sign with each node i, 7r(i) is the potential of node i. 
We define the reduced cost of arc of a set of node potentials 7r(i) 

cJ^:=Cij-n{i) + w{j). (48) 

The reduced costs in the residual network are defined in the same way as the 
costs, but with cjj instead of c^- . 
We have 

1) For any directed path P from k to J2{i,j)eP '^ij = J2{i,j)GP + d{k) - d{l). 

2) For any directed cycle W: E(i,j)ew <i = T.{i,j)ew Cij 

This means that negative cycles with respect to are also negative cycles with 
respect to c^j. 

Now we can formulate the reduced cost optimality condition: 

A feasible solution x* is an optimal solution of the min-cost flow problem 

<^ 

3 TT, a set of node potentials that satisfy the reduced cost optimality condition 
c7. >0 V(i,i) arc in G(x*). 

For the implication "<^" suppose that c^j > OV(i,j). One immediately real- 
izes that G(x*) contains no negative cycles since for each cycle W one has 
J2(i,j)ew '^ij — J2{i.j)ew '^ij ^ 0. For the other direction suppose that 

G(x*) contains no negative cycles. Denote with d{-) the shortest path distances 
from node 1 to all other nodes. Hence d{j) < d{i) + Cy S G(x*). Now 

define w = —d then cjj = Cij + d{i) — d{j) > 0. Note how closely connected the 
shortest path problem is to the min-cost flow problem. 

There is an intuitive economic interpretation of the reduced cost optimality 
condition. Suppose we interprete Cy as the cost of transporting one unit of a 
commodity from node i to node j through arc (i, j) and /Lt(i) = —n{i) as the 
cost of obtaining it at i. Then Cij + iJ,{i) is the cost of commdity at node j, if 
we obtain it at node i and transport it to node j via arc The inequality 

cjj > ix{j) < Cij + yu(z) says that the cost of commodity at node j is no 
more than obtaining it at i and sending it via — it could be smaller via 
other paths. 
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C.4 Successive shortest path algorithm 

With the concept of reduced costs we can now introduce the successive sliort- 
est patli algoritlim for solving the min-cost flow problem. The cycle canceling 
algorithm maintains feasil)ility of the solution at every step and attempts to 
achieve optimality. In contrast, the successive shortest path algorithm maintains 
optimality of the solution (c^ > 0) at every step and strives to attain feasibility 
(with respect to the mass balance constraints). 

A pseudoflow x: i?+ satisfies the capacity and non-negaivity constraints, 
but not necessarily the mass balance constraints. 
The imbalance of node i is defined as 

e{i} := b{i) + Yl - Yl ^iJ- (^9) 

OiaOeA} {j\Ui)eA} 

If e{i) > then we call e(z) the excess of node i, If e(i) < then we call it the 
deficit. E = {i\e{i) > 0} and E = {i\e{i) < 0} are the sets of excess and deficit 
nodes, respectively. Note that because of J2ieN ^(*) — J2ieN = we have 

E^ei^e(^) = -J2teD ^-i^)- 

Let the pseudoflow x satisfy the reduced cost optimality condition with re- 
spect to the node potential n and d{-) the shortest path distances from some 
node s to all the other nodes in the residual network G(x) with cjj as arc lengths. 
Therefore we have: 

Lemma 1: 

a) For the potential tt' = tt — li we have > 0, too. 

b) c^j = for all arcs (i, j) on shortest paths. 

To see a) note that d{j) < d{i) + cjj, thus c^- = Cij + {iT{i) — d{i)) — {Tr{j) — d{i)) = 
cjj + d{j) — d{i) > 0. For b) replace only the inequality by an equality. 

The following lemma is the basis of the siibsequent algorithm: Make the same 
assumption as in Lemma 1. Now send flow along a shortest path from some node 
s to some other node k to obtain a new pseudoflow x'. 

Lemma 2: 

x' also satisfies the reduced cost optimality conditions! 

For the proof take tt and ir' as in Lemma 1 and let P be the shortest path from 
node s to node k. Because of part b) of Lemma 1 it is ^ G P : c^- = 0. 
Therefore cj^ = —cjj = 0. Thus a flow augmentation on G P might add 
to the residual network, but cj^ = 0, which means that still the reduced 
cost optimality condition cj^ > is fulfilled. 

The strategy for an algorithm is now clear. By starting with a feasible solution 
that fulfills the reduced cost optimality condition one can iteratively send flow 
along the shortest paths from excess nodes to deflcit nodes to finally fulfill also 
the mass balance constraints. 
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algorithm successive shortest paths 
begin 

X := (G(x) =G) and TT := (cij = c^- > 0) 

e(i) := b{i) yi e N 

E := {i\e(i) > 0}, D := {i\e{i) < 0}. 

while ^^0 do 

begin 

select a node k G E and a node I £ D 

determine shortest path distance d{j) from some node s to 

all other nodes in G(x) w.r. to the reduced costs c^- 
let P denote a shortest path from node k to node I 
update TT := TT — d 

S := min {e(fc), — e(Z), min {r^ j) G P} } 
augment 6 units of flow along path P 
update X, (j(x), E, D and the reduced costs 



Note that in each iteration one excess is decreased by increasing flow ny at least 
one unit. Denoting with U the upper bound on the largest supply of any node 
one needs at most nil iterations, in each of which one has to solve a shortest path 
problem with non- negative arc lengths (so Dijkstra's algorithm is appropriate). 
This means that the above algorithm is polynomial if we know how U scales 
with m or n. 

C.5 Convex cost flows 

The cycle annealing algorithm as well as the successive shortest path algorithm 
solve the minimum cost flow problem for a linear cost function -^^^ Cij ■ Xij, 
where Cij represents the cost for sending one unit of flow along along the arc 
(i,j). This problem can be generalized to the following situation; 



subject to the mass balance constraints (46) and the capacity constraint (47). In 
addition we demand the flow variables Xij to be integer. The functions hij {xij ) 
can be any non-linear function, which has, however, to be convex, i.e. 



Va;, y , and ^ G [0, 1] hij{ex + (1 - e)y) < 9hij{x) + (1 - e)hij{y) (51) 



For this reason it is called convex cost flow problem. Without loss of generality 

we can assume that ft.y (0) = 0. Here the cost (for one unit) depends on the 
actual flow (since hij(xij) is a nonlinear function of the flow variable Xij): 



end 



end 



Minimize 




(50) 



(52) 
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Now Cij and cji are the costs for increasing and decreasing, respectively, the the 
flow variable Xij by one. 

After introducing these quantities it becomes straightforward to solve this 
problem with slight modifications of the algorithms we have already at hand. 
The first is again a negative cycle canceling algorithm: 

algorithm cycle canceling (convex costs) 
begin 

establish a feasible flow x 

calculate the costs c(x) as in eq. (52) 
while [G{x),c{x)] contains a negative cycle do 
begin 

use some algorithm to identify a negative cycle W 
augment one unit of flow in the cycle 
update G(x) and c{x) 
end 
end 

Note that since hij{x) is convex the cost for augmenting x by more than one 
unit increases the costs. This ensures that if we do not find any negative cycles, 
the flow is indeed optimal. 

This algorithm is, unfortunately non-polynomial in time, although it per- 
forms reasonably well on average. The successive shortest path algorithm dis- 
cussed in the last section can also be applied in the present context with a 
signiflcant gain in efliciency. For this algorithm it was essential that the reduced 
costs cjj with respect to some node potential tt maintains the reduced cost opti- 
mality condition c^j > upon flow augmentation along shortest paths. Now the 
question is, whether this still holds if with the change of the flow (caused by the 
augmentation) also the costs change. To show this we prove the folowing 
Lemma: Lot s be an excess node, d{-) shortest path distances w.r. to the reduced 
costs cjj from node s to all other nodes, w' = w + d, t a deficit node, W a shortest 
path from s to t, and x"*^ the flow that one obtains by augmenting x along W 
by one unit. Then: 

Cij > also for the modified arc costs along W . 

For the proof let Wij € {-|-1, — 1} for an arc € W with xlj"^ = Xij + Wij. 
Then the modified costs on this arc are 

because of the convexity of hij{x). From this follows for the modified reduced 
costs c^j* = c^j + 7r'(i) — 7r'(j) > Cij + 7r'(«) — 7r'(j), which proves the lemma. 

Thus we have the successive shortest path algorithm for the convex costs fiow 
problem: 
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algorithm successive shortest paths (convex costs) 
begin 

X := mill {i? (x) | x G Z^} and tt := 

e{i) := Hi) + T,{j\ui)eA} ^ji - E{j|Oi)eA} '^y Vi e 

while there is a node s with e(s) > do 

begin 

compute the reduced costs c'^(x) 

determine shortest path distance d{-) from s to 

all other nodes in (j(x) w.r. to the reduced costs cjj 
choose a node t with e(t) < 

augment x along shortest path from s to t by one unit 
IT = TT — d 
end 
end 

Note that we start with an optimal solution for -ff (x), i.e. we choose for each arc 
(z, j) the value of Xij in such a way that it is minimal. By this we guarantee that 
dj [xij ) > and thus that the reduced costs with tt = fulfill the optimality 
condition c^- > 0. The complexity of this algorithm is strictly polynomial in the 
physical example we discuss in section 8 [31]. 



Final remcirk 

For everybody who encounters one of the network flow problems mentioned 
above in his study of a physical (or any other) problem an important advice: 
Before reinventing the wheel, which means before wasting the time in hacking 
a subroutine that solves a particular network flow problem, one should consult 
the extremely valuable LEDA (library of efficient data types and algorithms) 
library, where many source codes of highly efficient combinatorial optimization 
algorithms can be found. All information, the manual [41] and the source codes 
can be found on the internet (this is public domain): 

http://www.mpi-sb.mpg.de/LEDA/leda.html 

Have fun! 
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